#############################################
## PROGRAM NAME: 202d_muni_regs.R          ##
## AUTHOR: MATT MLECZKO                    ##
## INPUTS:                                 ##
##         002_tract_to_pl.Rda             ##
##         002_tract_to_cs.Rda             ##
##         005_cc_out.Rda                  ##
##         201_af_muni_fin.Rda             ##
##                                         ##
## OUTPUTS:                                ##
##          202_muni_results_supp_1.qs -   ##
##          202_muni_results_supp_100.qs   ##
##                                         ##
## PURPOSE: Muni regressions               ##
##                                         ##
## LIST OF UPDATES:                        ##
#############################################

#log <- file(# USER DEFINED PATH AND FILE NAME HERE #)
#sink(log, append=TRUE)
#sink(log, append=TRUE, type="message")

set.seed(08542)

## load libraries ##

library(tidyverse)
library(foreign)
library(stringr)
library(tm)
library(gdata)
library(gsubfn)
library(haven)
library(readxl)
library(mice)
library(Rcpp)
library(basecamb)
library(qs)

## define paths ##
data_path <- # USER DEFINED PATH HERE #
pr_path <- # USER DEFINED PATH HERE #

## set working directory ##
setwd(data_path)

`%notin%` <- Negate(`%in%`)

## create a merge function that creates merge frequency as in Stata ##
## from user rwbuie at stackoverflow: https://stackoverflow.com/questions/30358401/is-there-a-way-to-create-statas-merge-indicator-variable-with-rs-merge ##
stata.merge <- function(x,y, name){
  x$df1 <- 1
  y$df2 <- 2
  df <- merge(x,y, by = name, all = TRUE)
  df$merge.variable <- rowSums(df[,c("df1", "df2")], na.rm=TRUE)
  df$df1 <- NULL
  df$df2<- NULL
  df
  #print(table(df$merge.variable))
  
  ## return the merged dataframe
  return(df)
}

## load analytic files ##

load("002_tract_to_pl.Rda")
load("002_tract_to_cs.Rda")

load("005_cc_out.Rda")
load("201_af_muni_fin.Rda")

af.muni.fin <- rename(af.muni.fin, 
                      rb_muniwide_2009 = rb_2009_muniwide,
                      rb_muniwide_2022 = rb_2022_muniwide)

## how extensive does the imputation need to be? ##

summary(af.muni.fin$zri_2006)
summary(af.muni.fin$sindex1_2006)
summary(af.muni.fin$sindex2_2006)
summary(af.muni.fin$sindex3_2006)
summary(af.muni.fin$sindex4_2006)
summary(af.muni.fin$sindex5_2006)
summary(af.muni.fin$sindex6_2006)
summary(af.muni.fin$sindex7_2006)

sapply(af.muni.fin, function(x) sum(is.na(x)))

## prep data for mi ##
af.muni.rd <- af.muni.fin %>%
  select(GEOID,
         pop_total_2009,
         pop_density_2009,
         per_18t34_2009,
         per_35t64_2009,
         per_65over_2009,
         median_age_2009,
         per_asian_2009,
         per_black_2009,
         per_latino_2009,
         per_white_2009,
         entropy_er_2009,
         per_child_2009,
         per_nomove_2009,
         per_immigrant_2009,
         hownrt_2009,
         median_gross_rent_2009,
         median_prop_value_2009,
         vacant_rt_2009,
         median_yr_str_2009,
         median_hhld_inc_2009,
         hhld_pov_rt_2009,
         unem_rt_2009,
         gini_2009,
         oo_hu_2009,
         per_bm_2009,
         sindex1_2006,
         sindex2_2006,
         sindex3_2006,
         sindex4_2006,
         sindex5_2006,
         sindex6_2006,
         sindex7_2006,
         sindex1_2022,
         sindex2_2022,
         sindex3_2022,
         sindex4_2022,
         sindex5_2022,
         sindex6_2022,
         sindex7_2022,
         zri_2006,
         zri_st_2006,
         zri_2022,
         zri_st_2022,
         zri_st_2006_median,
         zri_st_2022_median)

test <- af.muni.rd %>%
  select(-GEOID)

test.cor <- cor(test)

sapply(af.muni.rd, function(x) sum(is.na(x)))

## multiple imputation ## 
## source: https://library.virginia.edu/data/articles/getting-started-with-multiple-imputation-in-r ##

imp <- mice(af.muni.rd, maxit=0)

## extract predictorMatrix and methods of imputation ##

predM <- imp$predictorMatrix
meth <- imp$method
det(predM)

## leave these out of the predictor matrix ##
predM[, c("GEOID")] <- 0
predM[, c("sindex3_2006")] <- 0
predM[, c("sindex6_2006")] <- 0
predM[, c("sindex7_2006")] <- 0

## check determinant ##
det(predM)

## check the method ##
meth

## run the imputation ##

af.muni.mi1 <- mice(af.muni.rd, m = 10, 
                    predictorMatrix = predM, 
                    method = "cart", print =  FALSE)

## pool all the imputed data sets, along with original ## 

af.muni.complete <- complete(af.muni.mi1, "long", include = T)

summary(af.muni.complete)
sapply(af.muni.complete, function(x) sum(is.na(x)))


## PCA ## 

mi.pca <- function(impu){
  
  mi.pca.data <- af.muni.complete %>%
    filter(.imp == impu) %>%
    select(sindex1_2006, 
           sindex2_2006, 
           sindex3_2006, 
           sindex4_2006, 
           sindex5_2006, 
           sindex6_2006, 
           sindex7_2006)
  
  summary(mi.pca.data)
  
  mi.pca.res <- prcomp(mi.pca.data, center = T, scale=T)
  loadings <- mi.pca.res$rotation
  
  pca.data <- af.muni.complete %>%
    filter(.imp == impu) %>%
    mutate(zri_2006 = sindex1_2006*loadings[1,1] + 
             sindex2_2006*loadings[2,1] + 
             sindex3_2006*loadings[3,1] + 
             sindex4_2006*loadings[4,1] + 
             sindex5_2006*loadings[5,1] + 
             sindex6_2006*loadings[6,1] + 
             sindex7_2006*loadings[7,1])
  
  pca.data$zri_st_2006 <- (pca.data$zri_2006 - mean(pca.data$zri_2006, na.rm=T))/sd(pca.data$zri_2006,na.rm=T)
  
  return(pca.data)
  
}

mi.pca.out <- lapply(c(1,2,3,4,5,6,7,8,9,10), mi.pca)

## extract the pca output ## 

af.muni.out <- bind_rows(mi.pca.out)

summary(af.muni.out)

## get original data ##

af.muni.og <- af.muni.complete %>%
  filter(.imp == 0)

## bring on rest of covariates and outcomes ##

af.muni.m2 <- af.muni.fin %>%
  select(GEOID,
         pop_total_2022,
         pop_density_2022,
         per_18t34_2022,
         per_35t64_2022,
         per_65over_2022,
         median_age_2022,
         per_asian_2022,
         per_black_2022,
         per_latino_2022,
         per_white_2022,
         entropy_er_2022,
         per_child_2022,
         per_nomove_2022,
         per_immigrant_2022,
         hownrt_2022,
         median_gross_rent_2022,
         median_prop_value_2022,
         vacant_rt_2022,
         median_yr_str_2022,
         median_hhld_inc_2022,
         hhld_pov_rt_2022,
         unem_rt_2022,
         gini_2022,
         oo_hu_2022,
         per_bm_2022,
         starts_with("rb_"),
         starts_with("mpv_"),
         starts_with("mgr_"))

## merge ##
af.muni.reg.m <- stata.merge(af.muni.out,
                             af.muni.m2,
                             "GEOID")

## check merge ## 
table(af.muni.reg.m$merge.variable, useNA = "ifany")

## imputed data ##
af.muni.imp <- af.muni.reg.m %>%
  filter(merge.variable == 3) %>%
  select(-merge.variable) 

## now, original data ##

af.muni.m3 <- af.muni.fin %>%
  select(GEOID,
         pop_total_2022,
         pop_density_2022,
         per_18t34_2022,
         per_35t64_2022,
         per_65over_2022,
         median_age_2022,
         per_asian_2022,
         per_black_2022,
         per_latino_2022,
         per_white_2022,
         entropy_er_2022,
         per_child_2022,
         per_nomove_2022,
         per_immigrant_2022,
         hownrt_2022,
         median_gross_rent_2022,
         median_prop_value_2022,
         vacant_rt_2022,
         median_yr_str_2022,
         median_hhld_inc_2022,
         hhld_pov_rt_2022,
         unem_rt_2022,
         gini_2022,
         oo_hu_2022,
         per_bm_2022,
         starts_with("rb_"),
         starts_with("mpv_"),
         starts_with("mgr_"))

## merge ##
af.muni.og.fin.m <- stata.merge(af.muni.og,
                                af.muni.m3,
                                "GEOID")

## check merge ## 
table(af.muni.og.fin.m$merge.variable, useNA = "ifany")

## final reg file ##
af.muni.og.fin <- af.muni.og.fin.m %>%
  filter(merge.variable == 3) %>%
  select(-merge.variable)

## combine data ##

af.muni.cb <- rbind(af.muni.og.fin,
                    af.muni.imp)

af.muni.cb.fin <- af.muni.cb[order(af.muni.cb$.imp, af.muni.cb$.id),]

## final vars ##

af.muni.final <- af.muni.cb.fin %>%
  group_by(.imp) %>%
  mutate(zri_2006_fin = zri_st_2006 + zri_st_2006_median,
         zri_2022_fin = zri_st_2022 + zri_st_2022_median,
         zri_2006_fin_st = (zri_2006_fin - mean(zri_2006_fin, na.rm=T))/sd(zri_2006_fin,na.rm=T),
         zri_2022_fin_st = (zri_2022_fin - mean(zri_2022_fin, na.rm=T))/sd(zri_2022_fin,na.rm=T),
         treat_2006 = ifelse(zri_2006_fin_st > mean(zri_2006_fin_st, na.rm=T), 1,0),
         treat_2022 = ifelse(zri_2022_fin_st > mean(zri_2022_fin_st, na.rm=T), 1,0)) %>%
  ungroup() %>%
  mutate(zri.i = case_when(zri_2022 > zri_2006 + 0.5 ~ 1,
                           zri_2022 <= zri_2006 + 0.5 ~ 0),
         zri.d = case_when(zri_2022 < zri_2006 - 0.5 ~ 1,
                           zri_2022 >= zri_2006 - 0.5 ~ 0),
         zri.s = case_when(abs(zri_2022 - zri_2006) <= 0.5 ~ 1,
                           abs(zri_2022 - zri_2006) > 0.5 ~ 0),
         zri.id = abs(zri.i - zri.d),
         sindex1.change = sindex1_2022 - sindex1_2006,
         sindex2.change = sindex2_2022 - sindex2_2006,
         sindex3.change = sindex3_2022 - sindex3_2006,
         sindex4.change = sindex4_2022 - sindex4_2006,
         sindex5.change = sindex5_2022 - sindex5_2006,
         sindex6.change = sindex6_2022 - sindex6_2006,
         sindex7.change = sindex7_2022 - sindex7_2006,
         sindex1.i = case_when(sindex1_2022 > sindex1_2006 ~ 1,
                               sindex1_2022 <= sindex1_2006 ~ 0),
         sindex2.i = case_when(sindex2_2022 > sindex2_2006 ~ 1,
                               sindex2_2022 <= sindex2_2006 ~ 0),
         sindex3.i = case_when(sindex3_2022 > sindex3_2006 ~ 1,
                               sindex3_2022 <= sindex3_2006 ~ 0),
         sindex4.i = case_when(sindex4_2022 > sindex4_2006 ~ 1,
                               sindex4_2022 <= sindex4_2006 ~ 0),
         sindex5.i = case_when(sindex5_2022 > sindex5_2006 ~ 1,
                               sindex5_2022 <= sindex5_2006 ~ 0),
         sindex6.i = case_when(sindex6_2022 > sindex6_2006 ~ 1,
                               sindex6_2022 <= sindex6_2006 ~ 0),
         sindex7.i = case_when(sindex7_2022 > sindex7_2006 ~ 1,
                               sindex7_2022 <= sindex7_2006 ~ 0),
         sindex1.d = case_when(sindex1_2022 < sindex1_2006 ~ 1,
                               sindex1_2022 >= sindex1_2006 ~ 0),
         sindex2.d = case_when(sindex2_2022 < sindex2_2006 ~ 1,
                               sindex2_2022 >= sindex2_2006 ~ 0),
         sindex3.d = case_when(sindex3_2022 < sindex3_2006 ~ 1,
                               sindex3_2022 >= sindex3_2006 ~ 0),
         sindex4.d = case_when(sindex4_2022 < sindex4_2006 ~ 1,
                               sindex4_2022 >= sindex4_2006 ~ 0),
         sindex5.d = case_when(sindex5_2022 < sindex5_2006 ~ 1,
                               sindex5_2022 >= sindex5_2006 ~ 0),
         sindex6.d = case_when(sindex6_2022 < sindex6_2006 ~ 1,
                               sindex6_2022 >= sindex6_2006 ~ 0),
         sindex7.d = case_when(sindex7_2022 < sindex7_2006 ~ 1,
                               sindex7_2022 >= sindex7_2006 ~ 0),
         sindex1.s = case_when(sindex1_2022 == sindex1_2006 ~ 1,
                               sindex1_2022 != sindex1_2006 ~ 0),
         sindex2.s = case_when(sindex2_2022 == sindex2_2006 ~ 1,
                               sindex2_2022 != sindex2_2006 ~ 0),
         sindex3.s = case_when(sindex3_2022 == sindex3_2006 ~ 1,
                               sindex3_2022 != sindex3_2006 ~ 0),
         sindex4.s = case_when(sindex4_2022 == sindex4_2006 ~ 1,
                               sindex4_2022 != sindex4_2006 ~ 0),
         sindex5.s = case_when(sindex5_2022 == sindex5_2006 ~ 1,
                               sindex5_2022 != sindex5_2006 ~ 0),
         sindex6.s = case_when(sindex6_2022 == sindex6_2006 ~ 1,
                               sindex6_2022 != sindex6_2006 ~ 0),
         sindex7.s = case_when(sindex7_2022 == sindex7_2006 ~ 1,
                               sindex7_2022 != sindex7_2006 ~ 0))

summary(af.muni.final$zri.i)
summary(af.muni.final$zri.d)
summary(af.muni.final$zri.s)
summary(af.muni.final$zri.id)

summary(af.muni.final$sindex1.change)
summary(af.muni.final$sindex2.change)
summary(af.muni.final$sindex3.change)
summary(af.muni.final$sindex4.change)
summary(af.muni.final$sindex5.change)
summary(af.muni.final$sindex6.change)
summary(af.muni.final$sindex7.change)

summary(af.muni.final$sindex1.i)
summary(af.muni.final$sindex1.d)
summary(af.muni.final$sindex1.s)

summary(af.muni.final$sindex2.i)
summary(af.muni.final$sindex2.d)
summary(af.muni.final$sindex2.s)

summary(af.muni.final$sindex3.i)
summary(af.muni.final$sindex3.d)
summary(af.muni.final$sindex3.s)

summary(af.muni.final$sindex4.i)
summary(af.muni.final$sindex4.d)
summary(af.muni.final$sindex4.s)

summary(af.muni.final$sindex5.i)
summary(af.muni.final$sindex5.d)
summary(af.muni.final$sindex5.s)

summary(af.muni.final$sindex6.i)
summary(af.muni.final$sindex6.d)
summary(af.muni.final$sindex6.s)

summary(af.muni.final$sindex7.i)
summary(af.muni.final$sindex7.d)
summary(af.muni.final$sindex7.s)

## append original data ##
nzlu.panel.og <- af.muni.final %>%
  filter(.imp == 0) %>%
  mutate(zri_2006 = NA,
         zri.i = NA,
         zri.d = NA,
         zri.s = NA,
         zri.id = NA,
         sindex1.change = NA,
         sindex2.change = NA,
         sindex3.change = NA,
         sindex4.change = NA,
         sindex5.change = NA,
         sindex6.change = NA,
         sindex7.change = NA,
         sindex1.i = NA,
         sindex1.d = NA,
         sindex1.s = NA,
         sindex2.i = NA,
         sindex2.d = NA,
         sindex2.s = NA,
         sindex3.i = NA,
         sindex3.d = NA,
         sindex3.s = NA,
         sindex4.i = NA,
         sindex4.d = NA,
         sindex4.s = NA,
         sindex5.i = NA,
         sindex5.d = NA,
         sindex5.s = NA,
         sindex6.i = NA,
         sindex6.d = NA,
         sindex6.s = NA,
         sindex7.i = NA,
         sindex7.d = NA,
         sindex7.s = NA)

nzlu.panel.cb <- rbind(nzlu.panel.og,
                       af.muni.final)

nzlu.panel.mi2 <- as.mids(af.muni.final)
muni.nrows <- length(unique(af.muni.final$GEOID))

####################
## pooled results ## 
####################

## create no imp data ## 

af.muni.noimp.use <- af.muni.fin %>%
  mutate(zri_2006_fin_st = (zri_2006_cb - mean(zri_2006_cb, na.rm=T))/sd(zri_2006_cb,na.rm=T),
         zri_2022_fin_st = (zri_2022_cb - mean(zri_2022_cb, na.rm=T))/sd(zri_2022_cb,na.rm=T),
         treat_2006 = ifelse(zri_2006_fin_st > mean(zri_2006_fin_st, na.rm=T), 1,0),
         treat_2022 = ifelse(zri_2022_fin_st > mean(zri_2022_fin_st, na.rm=T), 1,0))

## create data for fixed effects regs ##

af.muni.fedata.t1 <- af.muni.final %>%
  select(GEOID,
         .imp,
         zri_2006_fin,
         ends_with("_2009")) %>%
  mutate(tp = "t1") %>%
  rename(zri_fin = zri_2006_fin) %>%
  rename_with(~str_remove(., '_2009'))

af.muni.fedata.t2 <- af.muni.final %>%
  select(GEOID,
         .imp,
         zri_2022_fin,
         pop_total_2022, 
         pop_density_2022,           
         per_18t34_2022,              
         per_35t64_2022,
         per_65over_2022,            
         median_age_2022,
         per_asian_2022,              
         per_black_2022,
         per_latino_2022,
         per_white_2022,
         entropy_er_2022,             
         per_child_2022, 
         per_nomove_2022,   
         per_immigrant_2022,           
         #oo_hu_2022,              
         hownrt_2022,
         median_gross_rent_2022,
         median_prop_value_2022,
         vacant_rt_2022,        
         median_yr_str_2022,
         median_hhld_inc_2022,  
         hhld_pov_rt_2022,          
         unem_rt_2022,
         gini_2022,
         oo_hu_2022,
         per_bm_2022,
         rb_npov_2022,
         rb_hpov_2022, 
         rb_hpov_cc_2022,
         rb_muniwide_2022,
         mgr_npov_2022,
         mgr_npov_rel_2022,
         mgr_hpov_2022, 
         mgr_hpov_rel_2022,
         mgr_hpov_cc_2022,
         mgr_hpov_cc_rel_2022,
         mpv_npov_2022,
         mpv_npov_rel_2022,
         mpv_hpov_2022, 
         mpv_hpov_rel_2022,
         mpv_hpov_cc_2022,
         mpv_hpov_cc_rel_2022) %>%
  mutate(tp = "t2") %>%
  rename(zri_fin = zri_2022_fin) %>%
  rename_with(~str_remove(., '_2022'))


af.muni.fedata <- rbind(af.muni.fedata.t1,
                        af.muni.fedata.t2)

muni.nrows.fe <- length(unique(af.muni.final$GEOID))*2

## create no imp data for fixed effects regs ##

af.muni.fedatanoimp.t1 <- af.muni.fin %>%
  select(GEOID,
         zri_2006_cb,
         ends_with("_2009")) %>%
  mutate(tp = "t1") %>%
  rename(zri_fin = zri_2006_cb) %>%
  rename_with(~str_remove(., '_2009')) 

af.muni.fedatanoimp.t2 <- af.muni.fin %>%
  select(GEOID,
         zri_2022_cb,
         pop_total_2022, 
         pop_density_2022,           
         land_area_2022,
         per_18t34_2022,              
         per_35t64_2022,
         per_65over_2022,            
         median_age_2022,
         per_asian_2022,              
         per_black_2022,
         per_latino_2022,
         per_white_2022,
         entropy_er_2022,             
         per_child_2022, 
         per_nomove_2022,   
         per_immigrant_2022,           
         oo_hu_2022,              
         hownrt_2022,
         median_gross_rent_2022,
         median_prop_value_2022,
         vacant_rt_2022,        
         median_yr_str_2022,
         median_hhld_inc_2022,  
         phi1_2022,
         phi2_2022,
         phi3_2022,
         phi4_2022,
         phi5_2022,
         phi6_2022,
         phi7_2022,
         phi8_2022,
         phi9_2022,
         phi10_2022,
         phi11_2022,
         phi12_2022,
         phi13_2022,
         phi14_2022,
         phi15_2022,
         phi16_2022,
         per_bm_2022,
         hhld_pov_rt_2022,          
         unem_rt_2022,
         gini_2022,
         rb_npov_2022,
         rb_hpov_2022, 
         rb_hpov_cc_2022,
         rb_muniwide_2022,
         mgr_npov_2022,
         mgr_npov_rel_2022,
         mgr_hpov_2022, 
         mgr_hpov_rel_2022,
         mgr_hpov_cc_2022,
         mgr_hpov_cc_rel_2022,
         mpv_npov_2022,
         mpv_npov_rel_2022,
         mpv_hpov_2022, 
         mpv_hpov_rel_2022,
         mpv_hpov_cc_2022,
         mpv_hpov_cc_rel_2022,) %>%
  mutate(tp = "t2") %>%
  rename(zri_fin = zri_2022_cb) %>%
  rename_with(~str_remove(., '_2022'))


af.muni.fedatanoimp.in <- rbind(af.muni.fedatanoimp.t1,
                                af.muni.fedatanoimp.t2) %>%
  mutate(zri_miss_flag = ifelse(is.na(zri_fin), 1, 0))

zri.miss.fe <- af.muni.fedatanoimp.in %>%
  filter(zri_miss_flag == 1)

af.muni.fedatanoimp <- af.muni.fedatanoimp.in %>%
  filter(GEOID %notin% zri.miss.fe$GEOID)

nrow(af.muni.fedatanoimp.in)
nrow(af.muni.fedatanoimp)

#####################
## begin bootstrap ##
#####################

source(paste0(pr_path,
              "msm_muni.R"))

source(paste0(pr_path,
              "msm_muni_noimp.R"))

## lists for storing output ##

muni.res1 <- list()
muni.res2 <- list()
muni.res3 <- list()
muni.res4 <- list()
muni.res5 <- list()
muni.res6 <- list()
muni.res7 <- list()
muni.res8 <- list()
muni.res9 <- list()

muni.res.noimp1 <- list()
muni.res.noimp2 <- list()
muni.res.noimp3 <- list()
muni.res.noimp4 <- list()
muni.res.noimp5 <- list()
muni.res.noimp6 <- list()
muni.res.noimp7 <- list()
muni.res.noimp8 <- list()
muni.res.noimp9 <- list()

tot_res <- function(i){
  
  ## muni level ##
  
  boot.munis <- data.frame(GEOID=sample(af.muni.final$GEOID,
                                        length(unique(af.muni.final$GEOID)),
                                        replace=T))
  
  boot.muni.data <- boot.munis %>%
    left_join(af.muni.final, "GEOID", relationship = "many-to-many") %>%
    arrange(.imp) %>%
    mutate(.id = rep(seq(1,muni.nrows, by=1),11)) %>%
    select(.imp,
           .id,
           GEOID,
           everything())
  
  boot.muni.noimp.data <- boot.munis %>%
    select(GEOID) %>%
    left_join(af.muni.noimp.use, "GEOID", relationship = "many-to-many") %>%
    ## drop missing ZRIs ##
    filter(!is.na(zri_2006_fin_st) & !is.na(zri_2022_fin_st))
  
  ## now, do the same for fixed effects data ## 
  
  boot.muni.fedata <- boot.munis %>%
    left_join(af.muni.fedata, "GEOID", relationship = "many-to-many") %>%
    arrange(.imp) %>%
    mutate(.id = rep(seq(1,muni.nrows.fe, by=1),11)) %>%
    select(.imp,
           .id,
           GEOID,
           everything())
  
  boot.muni.noimp.fedata <- boot.munis %>%
    select(GEOID) %>%
    left_join(af.muni.fedatanoimp, "GEOID", relationship = "many-to-many") %>%
    ## drop missing ZRIs ##
    filter(!is.na(zri_fin))
  
  print("OBS IMP")
  print(nrow(boot.muni.data))
  print(nrow(boot.muni.fedata))
  print("OBS NO IMP")
  print(nrow(boot.muni.noimp.data))
  print(nrow(boot.muni.noimp.fedata))
  
  af.muni.mids <- as.mids(boot.muni.data)
  af.muni.fe.mids <- as.mids(boot.muni.fedata)
  
  ############################
  ## Muni level regressions ##
  ############################
  
  muni.res1 <- msm_muni("mgr_npov_2009","mgr_npov_2022","mgr_npov", af.muni.mids, af.muni.fe.mids, 1)
  muni.res2 <- msm_muni("mgr_hpov_2009","mgr_hpov_2022","mgr_hpov", af.muni.mids, af.muni.fe.mids, 2)
  muni.res3 <- msm_muni("mgr_hpov_cc_2009","mgr_hpov_cc_2022", "mgr_hpov_cc", af.muni.mids, af.muni.fe.mids, 3)
  muni.res4 <- msm_muni("rb_npov_2009","rb_npov_2022","rb_npov", af.muni.mids, af.muni.fe.mids, 4)
  muni.res5 <- msm_muni("rb_hpov_2009","rb_hpov_2022","rb_hpov", af.muni.mids, af.muni.fe.mids, 5)
  muni.res6 <- msm_muni("rb_hpov_cc_2009","rb_hpov_cc_2022","rb_hpov_cc", af.muni.mids, af.muni.fe.mids, 6)
  muni.res7 <- msm_muni("mpv_npov_2009","mpv_npov_2022","mpv_npov", af.muni.mids, af.muni.fe.mids, 7)
  muni.res8 <- msm_muni("mpv_hpov_2009","mpv_hpov_2022","mpv_hpov", af.muni.mids, af.muni.fe.mids, 8)
  muni.res9 <- msm_muni("mpv_hpov_cc_2009","mpv_hpov_cc_2022","mpv_hpov_cc", af.muni.mids, af.muni.fe.mids, 9)
  
  muni.res.noimp1 <- msm_muni_noimp("mgr_npov_2009","mgr_npov_2022","mgr_npov", boot.muni.noimp.data, boot.muni.noimp.fedata, 1)
  muni.res.noimp2 <- msm_muni_noimp("mgr_hpov_2009","mgr_hpov_2022", "mgr_hpov",boot.muni.noimp.data, boot.muni.noimp.fedata, 2)
  muni.res.noimp3 <- msm_muni_noimp("mgr_hpov_cc_2009","mgr_hpov_cc_2022", "mgr_hpov_cc", boot.muni.noimp.data, boot.muni.noimp.fedata, 3)
  muni.res.noimp4 <- msm_muni_noimp("rb_npov_2009","rb_npov_2022","rb_npov", boot.muni.noimp.data, boot.muni.noimp.fedata, 4)
  muni.res.noimp5 <- msm_muni_noimp("rb_hpov_2009","rb_hpov_2022","rb_hpov", boot.muni.noimp.data, boot.muni.noimp.fedata, 5)
  muni.res.noimp6 <- msm_muni_noimp("rb_hpov_cc_2009","rb_hpov_cc_2022","rb_hpov_cc", boot.muni.noimp.data, boot.muni.noimp.fedata, 6)
  muni.res.noimp7 <- msm_muni_noimp("mpv_npov_2009","mpv_npov_2022","mpv_npov", boot.muni.noimp.data, boot.muni.noimp.fedata, 7)
  muni.res.noimp8 <- msm_muni_noimp("mpv_hpov_2009","mpv_hpov_2022","mpv_hpov", boot.muni.noimp.data, boot.muni.noimp.fedata, 8)
  muni.res.noimp9 <- msm_muni_noimp("mpv_hpov_cc_2009","mpv_hpov_cc_2022","mpv_hpov_cc", boot.muni.noimp.data, boot.muni.noimp.fedata, 9)
  
  all.res <- list(muni.res1,
                  muni.res2,
                  muni.res3,
                  muni.res4,
                  muni.res5,
                  muni.res6,
                  muni.res7,
                  muni.res8,
                  muni.res9,
                  muni.res.noimp1,
                  muni.res.noimp2,
                  muni.res.noimp3,
                  muni.res.noimp4,
                  muni.res.noimp5,
                  muni.res.noimp6,
                  muni.res.noimp7,
                  muni.res.noimp8,
                  muni.res.noimp9)
  
  filename <- paste0(data_path,
                     "202_muni_results_supp_", i, ".qs")
  
  qsave(all.res, filename)
  
  
}

## run the function ## 


lapply(1:100, tot_res)


## END OF PROGRAM ##

#sink()